function InternalForces(sols::Array{Float64, 2}, ie::Int64, step::Int64, point::Int64)
   sol = sols[step,:]
   num = gnum[ie,:]
   ecoord = gcoord[num,:]
   x1 = ecoord[1, 1]
   x3 = ecoord[3, 1]
   g = g_g[: , ie]
   a1 = -0.5*hweb - fbot; b1 = -0.5*hweb;
   a2 = -0.5*hweb; b2 = 0.5*hweb;
   a3 = 0.5*hweb; b3 = 0.5*hweb + ftop;
   D0 = D_matrixs(Es, μs)
   de = zeros(ndof)
   for k = 1:ndof
      if g[k] != 0
         de[k] = sol[g[k]]
      else
         de[k] = 0
      end
   end
   ξ = gaussdatax[point, 1]; x = 0.5*(1 - ξ)*x1 + 0.5*(1 + ξ)*x3;
   Mc = 0; Fnc = 0; Fc = 0; Ms = 0; Fns = 0; Fs = 0;

   for j = 1:ngpz
       zc = 0.5*Hc*gaussdataz[j, 1]
       wc = 0.5*Hc*gaussdataz[j, 2]
       nγc = N_γc(ξ, x1, x3)
       fnc = fem_dot(StressHistory[step, ie, :], Nstress(ξ, zc))
       mc = fnc*zc;
       Gc = beta_E(time0[step])*E28/(2*(1 + μc));
       fc = (5/6)*Gc*fem_dot(nγc, de);

       Mc = Mc + Bc*wc*mc
       Fnc = Fnc + Bc*wc*fnc
       Fc = Fc + Bc*wc*fc
    end

    frebar = Ar*(Er*strainc(ie, sol, ξ, zr) - fem_dot(StressHistory[step, ie, :], Nstress(ξ, zr)));
    mc = frebar*zr;
    Mc = Mc + mc;
    Fnc = Fnc + frebar;

    for j = 1:ngpz
       zs = 0.5*(b1 - a1)*gaussdataz[j, 1] + 0.5*(a1 + b1)
       ws = 0.5*(b1 - a1)*gaussdataz[j, 2]
       B = B_matrix_s(ξ, zs, x1, x3)
       stress = D0*B*de
       Ms = Ms + wbot*ws*stress[1]*zs
       Fns = Fns + wbot*ws*stress[1]
       Fs = Fs + wbot*ws*stress[3]
    end

   for j = 1:ngpz
       zs = 0.5*(b2 - a2)*gaussdataz[j, 1] + 0.5*(a2 + b2)
       ws = 0.5*(b2 - a2)*gaussdataz[j, 2]
       B = B_matrix_s(ξ, zs, x1, x3)
       stress = D0*B*de
       Ms = Ms + wweb*ws*stress[1]*zs
       Fns = Fns + wweb*ws*stress[1]
       Fs = Fs + wweb*ws*stress[3]
    end

   for j = 1:ngpz
      zs = 0.5*(b3 - a3)*gaussdataz[j, 1] + 0.5*(a3 + b3)
      ws = 0.5*(b3 - a3)*gaussdataz[j, 2]
      B = B_matrix_s(ξ, zs, x1, x3)
      stress = D0*B*de
      Ms = Ms + wtop*ws*stress[1]*zs
      Fns = Fns + wtop*ws*stress[1]
      Fs = Fs + wtop*ws*stress[3]
   end
   return [x, Mc, Fnc, Fc, Ms, Fns, Fs]
end

function post1(step::Int64)
    IntForceData = zeros((ngpx + 2)*nels, 7)
    # 存储高斯点处的内力数据
    for ie = 1:nels
        num = gnum[ie,:];
        ecoord = gcoord[num,:];
        x1 = ecoord[1, 1];
        x3 = ecoord[3, 1];
        for point = 1:ngpx
            IntForceData[1 + (ie - 1)*(ngpx + 2) + point, :] = InternalForces(sols, ie, step, point);
        end
        IntForceData[[1 + (ie - 1)*(ngpx + 2), ie*(ngpx + 2)], 1] = [x1, x3]
    end

    for ie = 1:nels
        left = lagrange(gaussdatax[: , 1], -1)
        right = lagrange(gaussdatax[: , 1], 1)
        d0 = IntForceData[(ie - 1)*(ngpx + 2) + 2 : (ie - 1)*(ngpx + 2) + 1 + ngpx, 2:end]
        IntForceData[ 1 + (ie - 1)*(ngpx + 2), 2:end] = d0'*left
        IntForceData[ ie*(ngpx + 2), 2:end] = d0'*right
    end

    return IntForceData
end


function stresses(ie::Int64, step::Int64)
   #计算单元g高斯点处的应力
   sol = sols[step,:]
   num = gnum[ie,:]
   ecoord = gcoord[num,:]
   x1 = ecoord[1, 1]
   x3 = ecoord[end, 1]
   g = g_g[:, ie]
   a1 = -0.5*hweb - fbot; b1 = -0.5*hweb;
   a2 = -0.5*hweb; b2 = 0.5*hweb;
   a3 = 0.5*hweb; b3 = 0.5*hweb + ftop;
   D0 = D_matrixs(Es, μs)
   de = zeros(ndof)
   for k = 1:ndof
      if g[k] != 0
         de[k] = sol[g[k]]
      else
         de[k] = 0
      end
   end
   sigmacx = zeros(ngpx, ngpz)
   taoc = zeros(ngpx, ngpz)

   sigmasx1 = zeros(ngpx, ngpz)
   sigmasz1 = zeros(ngpx, ngpz)
   taos1 = zeros(ngpx, ngpz)

   sigmasx2 = zeros(ngpx, ngpz)
   sigmasz2 = zeros(ngpx, ngpz)
   taos2 = zeros(ngpx, ngpz)

   sigmasx3 = zeros(ngpx, ngpz)
   sigmasz3 = zeros(ngpx, ngpz)
   taos3 = zeros(ngpx, ngpz)
for point = 1:ngpx
    ξ = gaussdatax[point, 1];
    x = 0.5*(1 - ξ)*x1 + 0.5*(1 + ξ)*x3;
   for j = 1:ngpz
      zc = 0.5* Hc*gaussdataz[j, 1]
      nγc = N_γc(ξ, x1, x3)
      Gc = beta_E(time0[step])*E28/(2*(1 + μc))
      sigmacx[point, j] = fem_dot(StressHistory[step, ie, :], Nstress(ξ, zc))
      taoc[point, j] = (5/6)*Gc*fem_dot(nγc, de);
   end
    for j = 1:ngpz
     zs = 0.5*(b1 - a1)*gaussdataz[j, 1] + 0.5*(a1 + b1);
     B = B_matrix_s(ξ, zs, x1, x3);
     stress = D0*B*de;
     sigmasx1[point, j] = stress[1];
     sigmasz1[point, j] = stress[2];
     taos1[point, j] = stress[3];
    end
    for j = 1:ngpz
     zs = 0.5*(b2 - a2)*gaussdataz[j, 1] + 0.5*(a2 + b2);
     B = B_matrix_s(ξ, zs, x1, x3);
     stress = D0*B*de;
     sigmasx2[point, j] = stress[1];
     sigmasz2[point, j] = stress[2];
     taos2[point, j] = stress[3];
    end
    for j = 1:ngpz
     zs = 0.5*(b3 - a3)*gaussdataz[j, 1] + 0.5*(a3 + b3);
     B = B_matrix_s(ξ, zs, x1, x3);
     stress = D0*B*de;
     sigmasx3[point, j] = stress[1];
     sigmasz3[point, j] = stress[2];
     taos3[point, j] = stress[3];
  end
end
   x = 0.5*(x3 - x1)*gaussdatax[: , 1] .+ 0.5*(x3 + x1);
   return (x, sigmacx, taoc, sigmasx1, sigmasz1, taos1, sigmasx2, sigmasz2, taos2, sigmasx3, sigmasz3, taos3)
end

function stress_for_originlab(step::Int64)
# 生成mathematica和originlab需要的等势图数据
    a1 = -0.5*hweb - fbot; b1 = -0.5*hweb;
    a2 = -0.5*hweb; b2 = 0.5*hweb;
    a3 = 0.5*hweb; b3 = 0.5*hweb + ftop;

    datacx = zeros(nels*ngpx*ngpz,3);
    datacxz = zeros(nels*ngpx*ngpz,3);

    datasx1 = zeros(nels*ngpx*ngpz,3);
    datasz1 = zeros(nels*ngpx*ngpz,3);
    datasxz1 = zeros(nels*ngpx*ngpz,3);

    datasx2 = zeros(nels*ngpx*ngpz,3);
    datasz2 = zeros(nels*ngpx*ngpz,3);
    datasxz2 = zeros(nels*ngpx*ngpz,3);

    datasx3 = zeros(nels*ngpx*ngpz,3);
    datasz3 = zeros(nels*ngpx*ngpz,3);
    datasxz3 = zeros(nels*ngpx*ngpz,3);

    count = 0;
    for iel = 1:nels
        (x, sigmacx, taoc, sigmasx1, sigmasz1, taos1, sigmasx2, sigmasz2, taos2, sigmasx3, sigmasz3, taos3) = stresses(iel, step);
        for sec = 1:ngpx
            for iz = 1:ngpz
                count = count + 1;

                datacx[count,:] = [x[sec], 0.5*Hc*gaussdataz[iz,1], sigmacx[sec, iz]]
                datacxz[count,:] = [x[sec], 0.5*Hc*gaussdataz[iz,1], taoc[sec, iz]]

                datasx1[count,:] = [x[sec], 0.5*(b1 - a1)*gaussdataz[iz, 1] + 0.5*(a1 + b1), sigmasx1[sec, iz]]
                datasz1[count,:] = [x[sec], 0.5*(b1 - a1)*gaussdataz[iz, 1] + 0.5*(a1 + b1), sigmasz1[sec, iz]]
                datasxz1[count,:] = [x[sec], 0.5*(b1 - a1)*gaussdataz[iz, 1] + 0.5*(a1 + b1), taos1[sec, iz]]

                datasx2[count,:] = [x[sec], 0.5*(b2 - a2)*gaussdataz[iz, 1] + 0.5*(a2 + b2), sigmasx2[sec, iz]]
                datasz2[count,:] = [x[sec], 0.5*(b2 - a2)*gaussdataz[iz, 1] + 0.5*(a2 + b2), sigmasz2[sec, iz]]
                datasxz2[count,:] = [x[sec], 0.5*(b2 - a2)*gaussdataz[iz, 1] + 0.5*(a2 + b2), taos2[sec, iz]]

                datasx3[count,:] = [x[sec], 0.5*(b3 - a3)*gaussdataz[iz, 1] + 0.5*(a3 + b3), sigmasx3[sec, iz]]
                datasz3[count,:] = [x[sec], 0.5*(b3 - a3)*gaussdataz[iz, 1] + 0.5*(a3 + b3), sigmasz3[sec, iz]]
                datasxz3[count,:] = [x[sec], 0.5*(b3 - a3)*gaussdataz[iz, 1] + 0.5*(a3 + b3), taos3[sec, iz]]
            end
        end
    end
    return  (datacx, datacxz, datasx1, datasz1, datasxz1, datasx2, datasz2, datasxz2, datasx3, datasz3, datasxz3)
end


function stress_for_matplotlib(step::Int64)
    # 生成matplotlib.pyplot需要的等势图数据
    a1 = -0.5*hweb - fbot; b1 = -0.5*hweb;
    a2 = -0.5*hweb; b2 = 0.5*hweb;
    a3 = 0.5*hweb; b3 = 0.5*hweb + ftop;

    datacx = zeros(nels*ngpx, ngpz);
    datacxz = zeros(nels*ngpx, ngpz);

    datasx2 = zeros(nels*ngpx, ngpz);
    datasz2 = zeros(nels*ngpx, ngpz);
    datasxz2 = zeros(nels*ngpx, ngpz);

    xp = zeros(nels*ngpx)

    for iel = 1:nels
        (x, sigmacx, taoc, sigmasx1, sigmasz1, taos1, sigmasx2, sigmasz2, taos2, sigmasx3, sigmasz3, taos3) = stresses(iel, step)
        for iz = 1:ngpz
            for sec = 1:ngpx
                col = (iel - 1)*ngpx + sec
                row = iz
                xp[col] = x[sec]
                datacx[col, row] = sigmacx[sec, iz]
                datacxz[col, row] = taoc[sec, iz]
                datasx2[col, row] = sigmasx2[sec, iz]
                datasz2[col, row] = sigmasz2[sec, iz]
                datasxz2[col, row] = taos2[sec, iz]
            end
        end
    end
    zpc = 0.5*Hc*gaussdataz[:,1];
    zps = 0.5*(b2 - a2)*gaussdataz[:,1] .+ 0.5*(b2 + a2)
    return  (xp, zpc, zps, datacx, datacxz, datasx2, datasz2, datasxz2)
end
